(library (grid-gen grid)

(export gr:new)

(import (guile);;(chezscheme)
        (rnrs)
        (grid-gen main)
        (grid-gen vars)
        (grid-gen code)
        (grid-gen math)
        (grid-gen symath)
        (grid-gen phys)
        (grid-gen HD)
        (grid-gen geom)
        (except (srfi srfi-1) map)
)

(define grid-vars '(rho p (u v w) (x y z)))

(define gr:new (lambda (dims geomf)
  (let* ((maxdim (apply max dims))
         (ndims (length dims))
         (dims1 (map (lambda (x) (+ (* 2 (comp-order)) x)) dims))
         (grd  (comp-array (PHYS:vect ndims) dims1))
         (grd1 (comp-array (PHYS:vect ndims) dims1))
         (tq  (comp-array (comp-float "tq") (cons (PHYS:qlen ndims) dims1)))
         (tcr (comp-float "cr"))
         (tau (comp-float "tau"))
         (err (comp-int "errflg" 0))
         (ptr (lambda (ofs)
           (map (lambda (x) `(+ ,(cur-ptr (car x)) ,(cadr x)))
                (zip (iota1 ndims) ofs)
           )
         ))
         (grd-v (lambda (nm)
           (lambda pa
             (if (not (null? pa))
                 (set-ind (vec-var grd nm) (car pa))
                 (vec-var grd nm)
             )
           )
         ))
         (zero-v (lambda p 0.0))
        )
    (grid:new1 dims
               (grd-v "rho")
               (grd-v "p")
               (grd-v "u")
               (if (> ndims 1)
                   (grd-v "v")
                   zero-v
               )
               (if (> ndims 2)
                   (grd-v "w")
                   zero-v
               )
               grd1
               tq
               tcr
               tau
               err
               ptr
               geomf
    )
  )
))

(define grid:new1 (lambda (dims v-rho v-p v-u v-v v-w
                                grd1 tq tcr tau0 err0 ptr geomf)
  (letrec*
    ((maxdim (apply max dims))
     (ndims (length dims))
     (zptr (map (lambda (x) 0) (iota1 ndims)))
     (i-co (lambda (i) (map (lambda (x) `(- ,x ,(comp-order))) i)))
     
     (dims2 (map (lambda (x) (+ x (* 2 (comp-order)))) dims))
     
     (geom (geom:new geomf))
     
     (err err0)
     (tau tau0)

     (t-rho (get-float-var))
     (t-p   (get-float-var))
     (t-u   (get-float-var))
     (t-v   (if (> ndims 1) (get-float-var) 0.0))
     (t-w   (if (> ndims 2) (get-float-var) 0.0))

     (add-d (lambda (l d o)
      (if (eq? d 0)
          (cons `(+ ,(car l) ,o)
                 (cdr l)
          )
          (cons (car l) (add-d (cdr l) (- d 1) o))
      )
     ))

     (bnd  (make-hash-table))
     (vbnd (list))
     
     (sbnd (lambda (i j f) (hash-set! bnd `(,i ,j) f)))

     (svbnd (lambda (typ pred f)
       (set! vbnd (append vbnd `((,typ ,pred ,f))))
;;       (format #t "SETVB: ~a\n" vbnd)
     ))

     (cofs (lambda (ar)
       (if (not (null? (cdr ar)))
           (cadr ar)
           zptr
       )
     ))
     
     (ipm (lambda (i d o) (replace-el i d (+ (list-ref i d) o))))
     
     (make-corr (lambda ()
       (beg
         (do-loop-in dims2 (comp-order)
           (apply beg
             (map (lambda (f) 
                   `(if ,(cons 'or corr-pred) ,(f self zptr) ())
                  )
                  corr
             )
           )
         )
         (raw "")
       )
     ))
     
     (mkbnd (lambda (d s cod)
       (do-loop `(,(comp-order))
         (do-loop-in-d dims2
                       (comp-order)
                       (filter (lambda (x) (not (= x d))) (iota1 ndims))
           (beg
             (set-ptr (list-ins (map cur-ptr (iota1 (dec ndims)))
                                d
                               `(+ ,(cur-ptr (dec ndims))
                                    (* ,s (+ ,(comp-order) ,(list-ref dims d)))
                                )
                      )
             )
             (apply beg cod)
           )
         )
       )
     ))

     (mk-flux (lambda (d gr i)
       (let* ((dF  (ADV:dF d gr i))
              (dF1 (ADV:dF1 d gr i))
              (Fl (gr 'Fl i d))
              (Fr (gr 'Fr i d))
              (ip1 (replace-el i d (+ (list-ref i d) 1)))
              
              (il (map (lambda (x)
                         (replace-el
                           i 
                           d
                           (+ (list-ref i d) (- x (dec (comp-order))))
                         )
                       )
                       (iota1 (* 2 (comp-order)))
                  )
              )

              (cndf (lambda (bnds ptrs)
                      (apply
                        append
                        (map (lambda (ptr)
                               (map (lambda (c)
                                      ((cadr c) gr ptr)
                                    )
                                    bnds
                               )
                             )
                             ptrs
                        )
                      )
                    )
              )

              (cnd (cndf vbnd il))
              (cnd-t (cons 'and cnd))
              (cnd-f (cons 'and (map (lambda (c) `(not ,c)) cnd)))
              
              (f-in-bnds (filter (lambda (b) (eq? (car b) 'free-in)) vbnd))
              (cnd-l (cons 'or (cndf f-in-bnds `(,i))))
              (cnd-r (cons 'or (cndf f-in-bnds `(,ip1))))

              (wall-bnds (filter (lambda (b) (eq? (car b) 'wall)) vbnd))
              (cnd-wl (cons 'or (cndf wall-bnds `(,i))))
              (cnd-wr (cons 'or (cndf wall-bnds `(,ip1))))
              
              (Fl0 (car (arr-elements Fl)))
              (Fr0 (car (arr-elements Fr)))
             )
         (map-ar (lambda (fl fr f2 f3)
                  `(if (or ,cnd-t ,cnd-f)
                       (- (* 0.5 (+ ,fl ,fr)) (+ ,f2 ,f3))
                       (if (or (and ,cnd-wl (not ,cnd-wr))
                               (and (not ,cnd-wl) ,cnd-wr)
                           )
                           0.0
                           (if (and (not ,cnd-l) ,cnd-r)
                               (if (>= ,Fl0 0.0) ,fl 0.0)
                               (if (and ,cnd-l (not ,cnd-r))
                                   (if (<= ,Fr0 0.0) ,fr 0.0)
                                   (- (* 0.5 (+ ,fl ,fr)) ,f2)
                               )
                           )
                       )
                   )
                 )
                 Fl Fr dF dF1
         )
       )
     ))
     
     (bkeys (make-hash-table))
     
     (self (lambda ar
             (letrec* 
               ((cmd (car ar))
                (spl-c (lambda (d)
                  (let* ((dm (cons ;; list of dims with d first
                               d
                               (filter (lambda (x) (not (= x d)))
                                       (iota1 ndims)
                               )
                             )
                         )
                         (cpi (map cur-ptr (iota1 ndims)))
                         (im  (replace-el zptr d (- (list-ref zptr d) 1)))
                         (xr (lambda (i) (self 'xcr i d)))
                         (xc (lambda (i) (self 'xc  i)))
                         (Qmc (apply (geom 'Qm) (xc zptr)))
                         (cf `(/ (* ,tau ,(+ (list-ref dims d) 0.0))
                                ,Qmc
                              )
                         )
                         (FF0 (lambda (i)
                           (mul2 cf
                                 (mul2 (apply (list-ref (geom 'gks) d) (xr i))
                                       (mk-flux d self i)
                                 )
                           )
                         ))
                         
                         (vr0 (self 'vr1-proj zptr (xr im) geom))

                         (vr0-c (PHYS:vrs:rpuvw t-rho t-p t-u t-v t-w))
                         
                         (q0 (PHYS:q vr0-c))
                         
                         (q1 (PHYS:q (PHYS:vrs-proj
                                       (PHYS:vrs-bkproj
                                         (PHYS:vrs:q
                                           (add q0 (FF0 im))
                                         )
                                         (xr im)
                                         geom
                                       )
                                       (xr zptr)
                                       geom
                                     )
                             )
                         )
                         (v2 (PHYS:vrs-bkproj
                               (PHYS:vrs:q (add q1 (mul2 -1.0 (FF0 zptr))))
                               (xr zptr)
                               geom
                             )
                         )
                         (v3 (if (= d (- ndims 1))
                                 (PHYS:srcs v2 (self 'xyz zptr) tau)
                                 v2
                             )
                         )
                        )
                    ((if (> ndims 1) omp-do beg)
                      (do-loop-in-d dims2 (comp-order) dm
                        (beg
                          (exp-sv
                            (set-ptr (list-ins (cdr cpi) d (cur-ptr 0)))
                            (omp-locals t-rho t-p t-u)
                            (if (> ndims 1) (omp-locals t-v) nop)
                            (if (> ndims 2) (omp-locals t-w) nop)
                            (set-val t-rho (vr0 'rho))
                            (set-val t-p (vr0 'p))
                            (set-val t-u (vr0 'u))
                            (if (> ndims 1) 
                                (set-val t-v (vr0 'v))
                                nop
                            )
                            (if (> ndims 2)
                                (set-val t-w (vr0 'w))
                                nop
                            )
                            
                            (set-val (self 'r1) (v3 'rho))
                            (set-val (self 'p1) (v3 'p))
                            (set-val (self 'u1) (v3 'u))
                            
                            (if (> ndims 1)
                                (set-val (self 'v1) (v3 'v))
                                nop
                            )
                            (if (> ndims 2)
                                (set-val (self 'w1) (v3 'w))
                                nop
                            )
                          )
                        )
                      )
                    )
                  )
                ))
               )
               
               (cond ((eq? cmd 'add-bnd)
                      (let ((ij (bnd2ij (cadr ar))))
                        (if ij
                            (sbnd (car ij) (cadr ij) (cddr ar))
                            (if (or (eq? (cadr ar) 'const)
                                    (eq? (cadr ar) 'free-in)
                                    (eq? (cadr ar) 'wall)
                                )
                                (svbnd (cadr ar) (caddr ar) (cadddr ar))
                                (#t (int-err "Invalid bound"))
                            )
                            
                        )
                      )
                     )

                     ((eq? cmd 'add-corr)
                      (set! corr
                            (append
                              corr
                              (list (cond ((eq? (cadr ar) 'smooth)
                                            (lambda (gr i)
                                              (cr-smooth)
                                            )
                                          )
                                          ((eq? (cadr ar) 'old)
                                            (lambda (gr i)
                                              (cr-old)
                                            )
                                          )
                                          ((eq? (cadr ar) 'err)
                                           (lambda (gr i)
                                             (cr-err)
                                           )
                                          )
                                          (#t (cadr ar))
                                    )
                              )
                            )
                      )
                     )

                     ((eq? cmd 'add-corr-pred)
                      (set! corr-pred
                            (cons
                              (cond ((eq? (cadr ar) 'tmax)
                                    `(>= ,(self 'p1) 
                                          (* ,PHYS:Rgas
                                             ,(self 'r1)
                                             ,(caddr ar)
                                          )
                                     )
                                    )
                                    ((eq? (cadr ar) 'rhomax)
                                    `(>= ,(self 'r1)
                                         ,(caddr ar)
                                     )
                                    )
                                    ((eq? (cadr ar) 'vmax)
                                    `(>= (+ (** ,(self 'u1) 2) 
                                            (** ,(self 'v1) 2)
                                            (** ,(self 'w1) 2)
                                         )
                                         (** ,(caddr ar) 2)
                                     )
                                    )
                                    (#t (cadr ar))
                              )
                              corr-pred
                            )
                      )
                     )
                     
                     ((eq? cmd 'advection)
                      (beg
                        (apply beg (map spl-c (iota1 ndims)))
                        (raw "")
                      )
                     )
               
                     ((eq? cmd 'c)  ((self 'var  (cofs ar)) 'c))


                     ((eq? cmd 'CFL)
     ;;;                 (format #t "CFL: ~a\n" ar)
                      (beg
                        (set-val tcr 0.0)
                        ((if (> ndims 1) omp-do beg)
                          (omp-reduce 'max tcr)
                          (do-loop dims2
                            (exp-sv
                              (set-val 
                                tcr
                                (let* ((vr (self 'var zptr))
                                       (cc (self 'c zptr))
                                       (vp ((geom 'proj) (math-cvector
                                                           (self 'u)
                                                           (self 'v)
                                                           (self 'w)
                                                         )
                                                         (self 'xc)
                                           )
                                       )
                                       (vc (map (lambda (x) 
                                                  `(+ (abs (* ,(car x)
                                                              ,(apply 
                                                                 (cadr x)
                                                                 (self 'xc)
                                                               )
                                                           )
                                                      ) 
                                                     ,cc
                                                   )
                                                )
                                                (zip (arr-elements vp)
                                                     (geom 'gks)
                                                )
                                           )
                                       )
                                       (vvl (map
                                              (lambda (x)
                                               `(* ,(list-ref vc x)
                                                   ,(+ 0.0 (list-ref dims x))
                                                )
                                              )
                                              (iota1 ndims)
                                            )
                                       )
                                     )
                                 `(max ,tcr
                                        (/ ,(cons 'max vvl)
                                           ,(apply (geom 'Qm) (self 'xc zptr))
                                        )
                                  )
                                )
                              )
                            )
                          )
                        )  
                        (set-val tau `(min ,tau (/ ,(cadr ar) ,tcr)))
                        (raw "")
                      )
                     )

                     ((eq? cmd 'correction)
                      (beg
                        (make-corr)
                        (raw "")
                      )
                     )

                     ((eq? cmd 'dims) dims)
                     ((eq? cmd 'dims-full) dims)

                     ((eq? cmd 'err) err)

                     ((eq? cmd 'F)
                      (PHYS:F (self 'var (cofs ar)))
                     )

                     ((eq? cmd 'F-proj)
                      (PHYS:F (self 'vr-proj (cofs ar) (caddr ar)))
                     )

                     ((eq? cmd 'Fl)   ;; (grid 'Fl ptr d) left F-proj for i-th right bound
                      ((PHYS:F (caddr ar)) (self 'var-l (cadr ar) (caddr ar)))
                     )

                     ((eq? cmd 'forall-grid-local)
                      (apply self
                        (append
                          `(forall-grid-local-d ,(iota1 ndims))
                           (cdr ar)
                        )
                      )
                     )

                     ((eq? cmd 'forall-grid-local-d)
                      (let ((dms (cadr ar))
                            (vrs (caddr ar))
                            (lam (cadddr ar))
                           )
                        (do-loop-in-d dims2 (comp-order) dms
                          (set-ptr (map cur-ptr dms))
                          (apply lam (map self vrs))
                        )
                      )
                     )

                     ((eq? cmd 'forall-grid)
                      (apply self (cons 'forall-grid-local (cdr ar)))
                     )
                     
                     ((eq? cmd 'forall-grid-d)
                      (apply self (cons 'forall-grid-local-d (cdr ar)))
                     )

                     ((eq? cmd 'Fr)   ;; (grid 'Fr ptr d) right F-proj for i-th right bound
                      ((PHYS:F (caddr ar)) (self 'var-r (cadr ar) (caddr ar)))
                     )

                     ((eq? cmd 'getbnd) (hash-ref bnd (cdr ar) #f))
                     
                     ((eq? cmd 'get-bnd-key)
                      (hash-ref bkeys (cadr ar) #f)
                     )
                     
                     ((eq? cmd 'get-vol-bnd)
                      vbnd
                     )

                     ((eq? cmd 'geom) geom)

                     ((eq? cmd 'init)
                      (beg
                        (do-loop dims2
                          ((cadr ar) self zptr)
                        )
                        (raw "")
                      )
                     )

                     ((eq? cmd 'init-grid)
                      nop
                     )
                     
                     ((eq? cmd 'load-grid-vals1)
                      (beg
                        (do-loop-in dims2 (comp-order)
                          (apply beg
                            (map (lambda (v1v2) 
                                   (apply set-val (map self v1v2))
                                 )
                                 (zip (self 'vrs1) (self 'vrs))
                            )
                          )
                        )
                        (set-val (self 'err) 0)
                        (raw "")
                      )
                     )

                     ((eq? cmd 'make-corr) make-corr)

                     ((eq? cmd 'mkbnd) (apply mkbnd (cdr ar)))

                     ((eq? cmd 'ndims) ndims)
                     
                     ((eq? cmd 'p) (v-p (ptr (cofs ar))))
                     ((eq? cmd 'p1)
                      (set-ind (vec-var grd1 "p") (ptr (cofs ar)))
                     )
                     ((eq? cmd 'pot)
                      (apply PHYS:Potential (self 'xyz (cofs ar)))
                     )
                     
                     ((eq? cmd 'ptr) (ptr (cofs ar)))
                     
                     ((eq? cmd 'q) (PHYS:q (self 'var (cofs ar))))

                     ((eq? cmd 'q-proj)
                      (PHYS:q (self 'vr-proj (cofs ar) (caddr ar)))
                     )

                     ((eq? cmd 'ql)   ;; (grid 'ql ptr d) left q-proj for i-th right bound
                      (PHYS:q (self 'var-l (cadr ar) (caddr ar)))
                     )

                     ((eq? cmd 'qr)   ;; (grid 'qr ptr d) right q-proj for i-th right bound
                      (PHYS:q (self 'var-r (cadr ar) (caddr ar)))
                     )

                     ((eq? cmd 'r) (v-rho (ptr (cofs ar))))
                     ((eq? cmd 'r1)
                      (set-ind (vec-var grd1 "rho") (ptr (cofs ar)))
                     )

                     ((eq? cmd 'save-all)
                      (beg
                        (print "Saving to file: " (cadr ar))
                        (to-file (cadr ar)
                          (do-loop-in dims (comp-order)
                            (apply fprint
                              (append
                                (list-head (self 'xyz) ndims)
                                (list (self 'r) (self 'p))
                                (list-head (self 'vels) ndims)
                              )
                            )
                          )
                        )
                      )
                     )
                     
                     ((eq? cmd 'self) self)
                     
                     ((eq? cmd 'set-bnd-key)
                      (hash-set! bkeys (cadr ar) (caddr ar))
                     )
                     
                     ((eq? cmd 'set-err) (set! err (cadr ar)))

                     ((eq? cmd 'set-flux) (set! mk-flux (cadr ar)))

                     ((eq? cmd 'set-geom) (set! geom (geom:new (cadr ar))))

                     ((eq? cmd 'set-make-corr) (set! make-corr (cadr ar)))

                     ((eq? cmd 'set-tau) (set! tau (cadr ar)))
                     
                     ((eq? cmd 'subgrids) `(,self))

                     ((eq? cmd 'tau) tau)
                     
                     ((eq? cmd 'temp)
                      ((self 'var (cofs ar)) 'temp)
                     )
                     
                     ((eq? cmd 'u) (v-u (ptr (cofs ar))))
                     ((eq? cmd 'u1)
                      (set-ind (vec-var grd1 "u") (ptr (cofs ar)))
                     )

                     ((eq? cmd 'update-bnd)
                      (beg
                        (apply beg
                          (hash-map->list
                            (lambda (ds cd)
                              (if (and (= (length cd) 5)
                                       (eq? (car cd) 'grid)
                                  )
                                  (let* ((g2 (list-ref cd 1))
                                         (ij (bnd2ij (list-ref cd 2)))
                                         (vr (list-ref cd 3))
                                         (mp (list-ref cd 4))
                                         (pt (lambda (s)
                                               (map
                                                 (lambda (d)
                                                  `(- ,(cur-ptr d)
                                                      ,(if (= d (car ij))
                                                           (if (= s 0)
                                                               0
                                                               (+ (comp-order)
                                                                  (list-ref
                                                                    dims
                                                                    (car ij)
                                                                  )
                                                               )
                                                           )
                                                           (comp-order)
                                                       )
                                                   )
                                                 )
                                                 (iota1 ndims)
                                               )
                                             )
                                         )
                                        )
                                    (beg
                                      (g2 'mkbnd (car ij) (cadr ij)
                                        (map (lambda (v)
                                               (set-val
                                                 (set-ind (cadr v) 
                                                          (pt (cadr ij))
                                                 )
                                                 (g2 (car v)
                                                     (replace-el
                                                       zptr
                                                       (car ij)
                                                       (if (= 0 (cadr ij))
                                                           (comp-order)
                                                           (* -1 (comp-order))
                                                       )
                                                     )
                                                 )
                                               )
                                             )
                                             (zip (g2 'vrs) vr)
                                        )
                                      )
                                      (self 'mkbnd (car ds) (cadr ds)
                                        (map (lambda (v)
                                               (set-val
                                                 (set-ind
                                                   (del-ind (self (car v)))
                                                   (map
                                                     (lambda (d)
                                                       (cur-ptr
                                                         (let 
                                                           ((a (assoc d mp)))
                                                           (if a (cadr a) d)
                                                         )
                                                       )
                                                     )
                                                     (iota1 ndims)
                                                   )
                                                 )
                                                 (set-ind (cadr v)
                                                          (pt (- 1 (cadr ij)))
                                                 )
                                               )
                                             )
                                             (zip (self 'vrs) vr)
                                        )
                                      )
                                    )
                                  )
                                  (self 'mkbnd (car ds) (cadr ds)
                                        (cd self zptr)
                                  )
                              )
                            )
                            bnd
                          )
                        )
                        (self 'volbnd)
                        (raw "")
                      )
                     )

                     ((eq? cmd 'update-grid-vals)
                      (beg
                       `(if (= ,(self 'err) 0)
                           ,(do-loop-in dims2 (comp-order)
                              (let ((cod (apply beg
                                                (map (lambda (v1v2)
                                                       (apply set-val 
                                                              (map self v1v2)
                                                       )
                                                     )
                                                     (zip (self 'vrs)
                                                          (self 'vrs1)
                                                     )
                                                )
                                         )
                                    )
                                    (volb (cons 'or
                                                 (map (lambda (c)
                                                        (c self zptr)
                                                      )
                                                      (map cadr vbnd)
                                                 )
                                          )
                                    )
                                   )
                                (if (null? vbnd)
                                    cod
                                   `(if  (not ,volb)
                                        ,cod
                                         ()
                                    )
                                )
                              )
                            )     
                            ()
                        )
                      )
                     )

                     ((eq? cmd 'v) (v-v (ptr (cofs ar))))
                     ((eq? cmd 'v1) 
                      (if (> ndims 1)
                          (set-ind (vec-var grd1 "v") (ptr (cofs ar)))
                          0.0
                      )
                     )

                     ((eq? cmd 'var)
                      (PHYS:vrs:rpuvw  (self 'r (cofs ar))
                                       (self 'p (cofs ar))
                                       (self 'u (cofs ar))
                                       (self 'v (cofs ar))
                                       (self 'w (cofs ar))
                      )
                     )

                     ((eq? cmd 'var1)
                      (PHYS:vrs:rpuvw  (self 'r1  (cofs ar))
                                       (self 'p1  (cofs ar))
                                       (self 'u1  (cofs ar))
                                       (self 'v1  (cofs ar))
                                       (self 'w1  (cofs ar))
                      )
                     )

                     ((eq? cmd 'var-l)   ;; (grid 'var-l ptr d) left v-proj for i-th right bound
                      (self 'vr-proj
                             (cofs ar)
                             (self 'xcr (cofs ar) (caddr ar))
                      )
                     )

                     ((eq? cmd 'var-r)   ;; (grid 'var-r ptr d) right v-proj for i-th right bound
                      (self 'vr-proj
                             (ipm (cofs ar) (caddr ar) 1)
                             (self 'xcr (cofs ar) (caddr ar))
                      )
                     )
                     
                     ((eq? cmd 'vels)
                      (map (lambda (v) (self v (cofs ar)))
                          '(u v w)
                      )
                     )
                     
                     ((eq? cmd 'vlist) (list (self 'u) (self 'v) (self 'w)))
                     
                     ((eq? cmd 'vr-proj)
                      (PHYS:vrs-proj (self 'var (cofs ar)) (caddr ar) geom)
                     )

                     ((eq? cmd 'vr1-proj)
                      (PHYS:vrs-proj (self 'var1 (cofs ar)) (caddr ar) geom)
                     )
                     
                     
                     ((eq? cmd 'vrs) (list-head '(r p u v w) (+ 2 ndims)))

                     ((eq? cmd 'vrs1) (list-head '(r1 p1 u1 v1 w1) (+ 2 ndims)))
  
                     
                     ((eq? cmd 'volbnd)
                      (let ((bnds (filter (pip not null? caddr) vbnd)))
                        (if (not (null? bnds))
                            (do-loop dims2
                              (apply beg
                                (map (lambda (bnd)
                                       `(if ,((cadr bnd) self zptr)
                                            ,((caddr bnd) self zptr)
                                             ()
                                        )
                                     )
                                     bnds
                                )
                              )
                            )
                        )
                      )
                     )
  
                     ((eq? cmd 'w) (v-w (ptr (cofs ar))))
                     ((eq? cmd 'w1)
                      (if (> ndims 2)
                          (set-ind (vec-var grd1 "w") (ptr (cofs ar)))
                          0.0
                      )
                     )

                     ((eq? cmd 'x) 
                      (car (apply self (cons 'xyz (cdr ar))))
                     )

                     ((eq? cmd 'xc)
                      (map
                        (lambda (x)
                          `(/ (- ,(cadr x) 1.0)
                              (- ,(car x) 1.0)
                           )
                        )
                        (zip dims (ptr (i-co (cofs ar))))
                      )
                     )

                     ((eq? cmd 'xcr) ;; Invoke: (grid 'xcr ptr dim)
                      (map
                        (lambda (x)
                          `(/ (- ,(cadr x)
                                 ,(if (= (car x) (caddr ar)) 0.5 1.0)
                              )
                              (- ,(list-ref dims (car x)) 1.0)
                           )
                        )
                        (zip (iota1 ndims) (ptr (i-co (cofs ar))))
                      )
                     )

                     ((eq? cmd 'xyz)
                      (apply (geom 'f) (apply self (cons 'xc (cdr ar))))
                     )

                     ((eq? cmd 'y)
                      (cadr (apply self (cons 'xyz (cdr ar))))
                     )
                     ((eq? cmd 'z)
                      (caddr (apply self (cons 'xyz (cdr ar))))
                     )
                     
                    (#t (int-err (format #f "Invalid grid query '~a'" cmd)))
               )
             )
           )
     )
     (cr-smooth
       (lambda ()
         (letrec* ((r-mid (lambda (vl d o)
                     (if (< d (self 'ndims))
                         (if (<= o (comp-order))
                             (let* ((im1 (ipm zptr d -1))
                                    (ip1 (ipm zptr d  1))
                                   )
                              `(+ (* (/ 1.0
                                       ,(+ 1.0 (* 2.0
                                                  (comp-order)
                                                  (self 'ndims)
                                               )
                                        )
                                      )
                                      (+ ,(self vl im1) ,(self vl ip1))
                                  )
                                 ,(r-mid vl d (+ o 1))
                               )
                             )
                            (r-mid vl (+ d 1) 1)
                        )
                       `(* (/ 1.0 ,(+ 1 (* 2.0 (self 'ndims)))) ,(self vl zptr))
                     )
                   ))
                  )
           (beg
             (set-val (self 'r1) (r-mid 'r 0 1))
             (set-val (self 'p1) (r-mid 'p 0 1))
             (set-val (self 'u1) (r-mid 'u 0 1))
             (if (> (self 'ndims) 1)
               (set-val (self 'v1) (r-mid 'v 0 1))
               nop
             )
             (if (> (self 'ndims) 2)
               (set-val (self 'w1) (r-mid 'w 0 1))
               nop
             )
           )
         )
       )
     )
     
     (cr-old
       (lambda ()
         (beg
           (set-val (self 'r1) (self 'r))
           (set-val (self 'p1) (self 'p))
           (set-val (self 'u1) (self 'u))
           (if (> (self 'ndims) 1)
               (set-val (self 'v1) (self 'v))
               nop
           )
           (if (> (self 'ndims) 2)
               (set-val (self 'w1) (self 'w))
               nop
           )
         )
       )
     )
     
     (cr-err (lambda () (set-val (self 'err) 1)))

     (corr (list))

     (corr-pred (append `((isnan ,(self 'r1))
                          (isnan ,(self 'p1))
                          (<= ,(self 'r1) 0.0)
                          (<= ,(self 'p1) 0.0)
                          (isnan ,(self 'u1))
                         )
                         (if (> ndims 1)
                            `((isnan ,(self 'v1)))
                            '()
                         )
                         (if (> ndims 2)
                            `((isnan ,(self 'w1)))
                            '()
                         )
                )
     )
    )
    self
  )
))

(define gr:norm (lambda (gr cmd eqn)
  (lambda args
    (if (and (not (null? args)) (member (car args) cmd))
        (eqn (apply gr args))
        (apply gr args)
    )
  )
))



(define gr:join (lambda (g1 b1 vr1 g2 b2 vr2)
  (if (not (equal? (g1 'vrs) (g2 'vrs)))
      (int-err "Cannot join grids: var list mismatch!")
  )
  (let* ((ij1 (bnd2ij b1))
         (ij2 (bnd2ij b2))

         (dims1 (replace-el (g1 'dims) (car ij1) (comp-order)))
         (dims2 (replace-el (g2 'dims) (car ij2) (comp-order)))
         
         ;; FIXME: take axis maping into account!
         
         (chk (if (not (equal? dims1 dims2))
                  (int-err "Cannot join grids: dims mismatch!")
              )
         )
         
         (vrl (map (lambda (v) 
                     (comp-array (comp-float (string-append "tbnd")) dims1)
                   )
                   (g1 'vrs)
              )
         )
         
;;         (bk (lambda () #t))
         (bk (gensym))
         
         (v2i (lambda (vrs)
           (map (lambda (v) (list-index (lameq v) '(x y z))) vrs)
         ))
         (vi1 (v2i vr1))
         (vi2 (v2i vr2))
         (mp12 (zip vi1 vi2))
         (mp21 (zip vi2 vi1))
        )

    (g1 'set-bnd-key ij1 bk)
    (g2 'set-bnd-key ij2 bk)

    (g1 'add-bnd b1 'grid g2 b2 vrl mp12)
    (g2 'add-bnd b2 'grid g1 b1 vrl mp21)
  )
))

)

(set! grid:cart gr:new)
(set! grid:norm gr:norm)
(set! grid:join gr:join)
